//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// --------------------------------------------------------------
//   GEANT 4 - Underground Dark Matter Detector Advanced Example
//
//      For information related to this code contact: Alex Howard
//      e-mail: alexander.howard@cern.ch
// --------------------------------------------------------------
// Comments
//
//                  Underground Advanced
//               by A. Howard and H. Araujo 
//                    (27th November 2001)
//
// DetectorConstruction program
// --------------------------------------------------------------

#include "DMXDetectorConstruction.hh"
#include "DMXDetectorMessenger.hh"

#include "DMXScintSD.hh"
#include "DMXPmtSD.hh"

#include "G4Material.hh"
#include "G4MaterialTable.hh"
#include "G4Element.hh"
#include "G4Isotope.hh"
#include "G4UnitsTable.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Torus.hh"
#include "G4Sphere.hh"
#include "G4UnionSolid.hh"
#include "G4SubtractionSolid.hh"

#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"

//#include "MyDeviceParameterisation.hh"
//#include "G4PVParameterised.hh"

#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"
#include "G4Transform3D.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4OpBoundaryProcess.hh"

#include "G4FieldManager.hh"
#include "G4UniformElectricField.hh"
#include "G4TransportationManager.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4EqMagElectricField.hh"
#include "G4ClassicalRK4.hh"
#include "G4ChordFinder.hh"

#include "G4SDManager.hh"

#include "G4VisAttributes.hh"
#include "G4Colour.hh"

#include "G4UserLimits.hh"

#include "G4RunManager.hh"

#include <TROOT.h>
#include <TH1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TRandom.h>
#include <TNtuple.h>


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
DMXDetectorConstruction::DMXDetectorConstruction()  
{
  // create commands for interactive definition of time cuts:
  detectorMessenger = new DMXDetectorMessenger(this);

  theUserLimitsForRoom     = 0; 
  theUserLimitsForWorld    = 0; 
  theUserLimitsForDetector = 0; 
  // default time cut = infinite
  //  - note also number of steps cut in stepping action = MaxNoSteps
  theMaxTimeCuts      = DBL_MAX;
  theMaxStepSize      = DBL_MAX;
  theDetectorStepSize = DBL_MAX;
  theWorldStepSize    = 1000.0*mm;
  theRoomTimeCut      = 1000. * nanosecond;
  theMinEkine         = 250.0*eV; // minimum kinetic energy required in volume
  theRoomMinEkine     = 250.0*eV; // minimum kinetic energy required in volume

}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
DMXDetectorConstruction::~DMXDetectorConstruction() 
{
  delete theUserLimitsForRoom;
  delete theUserLimitsForWorld;
  delete theUserLimitsForDetector;
  delete detectorMessenger;
}



//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void DMXDetectorConstruction::DefineMaterials() 
{

#include "DMXDetectorMaterial.icc"

}

// Get Argon Scintillation Spectrum

G4double DMXDetectorConstruction::GetArScintillationSpectrum(G4double kk)
{
  G4double waveL;
  waveL =exp(-0.5*((kk-128.0)/(2.929))*((kk-128.0)/(2.929)));      //Gaussian Form
  return waveL;
}

//Calculates the dielectric constant of GAr from the Bideau-Sellmeier formula.
//See : A. Bideau-Mehu et al., "Measurement of refractive indices of Ne, Ar,
//Kr and Xe ...", J. Quant. Spectrosc. Radiat. Transfer, Vol. 25 (1981), 395

G4double DMXDetectorConstruction::GetLArEpsilon(G4double lambda)
{
  G4double epsilon;
  if (lambda < (lowLambda * nanometer)) return 1.0e4; // lambda MUST be > 110.0 nm
  epsilon = lambda / micrometer; // switch to micrometers 
  epsilon = 1.0 / (epsilon * epsilon); // 1 / (lambda)^2 
  epsilon = 1.2055e-2 * ( 0.2075 / (91.012 - epsilon) +
                          0.0415 / (87.892 - epsilon) +
                          4.3330 / (214.02 - epsilon) );
  epsilon *= (8./12.); // Bideau-Sellmeier -> Clausius-Mossotti
  epsilon *= (LArRho / GArRho); // density correction (Ar gas -> LAr liquid)
  if (epsilon < 0.0 || epsilon > 0.999999) return 4.0e6;
  epsilon = (1.0 + 2.0 * epsilon) / (1.0 - epsilon); // solve Clausius-Mossotti
  return epsilon;
}

G4double DMXDetectorConstruction::GetLArRefIndex(G4double lambda)
{
  return ( sqrt(GetLArEpsilon(lambda)) ); // square root of dielectric constant 
}

//Calculates the Rayleigh scattering length using equations given in 
//G. M. Seidel at al., "Rayleigh scattering in rare-gas liquids",
//arXiv:hep-ex/0111054 v2 22 Apr 2002

G4double DMXDetectorConstruction::GetLArRayLength(G4double lambda,G4double temp)
{
  G4double h;
  h = GetLArEpsilon(lambda);
  if (h < 1.00000001) h = 1.00000001; // just a precaution
  h = (h - 1.0) * (h + 2.0); // the "dielectric constant" dependance
  h *= h; // take the square 
  h *= LArKT * temp * k; // compressibility * temp * Boltzmann constant 
  h /= lambda * lambda * lambda * lambda; // (lambda)^4 
  h *= 9.18704494231105429; // (2 * Pi / 3)^3 
  if ( h < (1.0 / (10.0 * km)) ) h = 1.0 / (10.0 * km); // just a precaution 
  if ( h > (1.0 / (0.1 * nanometer)) ) h = 1.0 / (0.1 * nanometer); // just a precaution
  return ( 1.0 / h );
}

G4double DMXDetectorConstruction::GetLArAbsLength(G4double lambda)
{
  const G4double d = 5.8; // cm
  const G4double a = 0.234; // +/- 0.014 nm-1 
  const G4double b = 113.02; // +/- 0.65 nm
  G4double T =  1.-exp(-a*(lambda-b));
  return d/log(1./T);
}

/*
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void DMXDetectorConstruction::DefineField() {

#include "DMXDetectorField.icc"

}
*/


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4VPhysicalVolume* DMXDetectorConstruction::Construct() {

  DefineMaterials();

  // DefineField();

  // make colours
  G4Colour  white   (1.0, 1.0, 1.0) ;
  G4Colour  grey    (0.5, 0.5, 0.5) ;
  G4Colour  lgrey   (.85, .85, .85) ;
  G4Colour  red     (1.0, 0.0, 0.0) ;
  G4Colour  blue    (0.0, 0.0, 1.0) ;
  G4Colour  cyan    (0.0, 1.0, 1.0) ;
  G4Colour  magenta (1.0, 0.0, 1.0) ; 
  G4Colour  purple  (.75, 0.0, 1.0) ; 
  G4Colour  yellow  (1.0, 1.0, 0.0) ;
  G4Colour  orange  (.75, .55, 0.0) ;
  G4Colour  lblue   (0.0, 0.0, .75) ;
  G4Colour  lgreen  (0.0, .75, 0.0) ;
  G4Colour  green   (0.0, 1.0, 0.0) ;
  G4Colour  brown   (0.7, 0.4, 0.1) ;
  //  un-used colours:
  //  G4Colour  black   (0.0, 0.0, 0.0) ;

  //----------------------------Air Property--------------------------------------//

  G4double Air_PP[3]    = { 9.286*eV, 9.686*eV, 10.186*eV };
  G4double Air_RIND[3]  = { 1.0000, 1.0000, 1.0000 };
  //  G4double Air_RIND[3]  = { 0.0000, 1.0000, 0.0000 };
  G4double Air_ABSL[3]  = { 1000*m, 1000*m, 1000*m}; //atten length

  G4MaterialPropertiesTable *Air_mt = new G4MaterialPropertiesTable();
  Air_mt->AddProperty("RINDEX",        Air_PP, Air_RIND,  3);
  Air_mt->AddProperty("ABSLENGTH",     Air_PP, Air_ABSL,  3);
  lab_mat->SetMaterialPropertiesTable(Air_mt);

  //--------------------------Gas Argon Property-------------------------------//

  // G4double GAr_RINDEX = 1.0004;
  // G4double GAr_ABSORP = 1000000*m;

  // G4double GAr_PP[20]    = {   2.951*eV,   2.952*eV,   2.953*eV, 8.98626*eV, 8.98627*eV, 9.08627*eV, 9.18627*eV, 9.28627*eV, 9.38627*eV, 9.48627*eV, 9.58627*eV, 9.68627*eV, 9.78627*eV, 9.88627*eV, 9.98627*eV, 10.0863*eV, 10.1863*eV, 10.2863*eV,  10.3863*eV,  10.3864*eV };
  // //  G4double GAr_PP[6]    = { 1.9*eV, 2.0*eV, 2.1*eV, 4.0*eV, 6.0*eV, 8.0*eV };
  // //  G4double GAr_PP[6]    = { 9.286*eV, 9.686*eV, 10.186*eV };
  // G4double GAr_SCINT[20] = {         0.,         0.,         0.,         0., 0.00339483,  0.0166968,  0.0636808,   0.180736,   0.394466,    0.67708,   0.918868,        1.0,   0.880227,   0.650421,   0.398116,   0.207939,  0.0927144,  0.0355237,   0.0115269,          0. };
  // //  G4double GAr_SCINT[6] = {  0.000,  0.000,  0.000,  0.0,  0.000,   0.0 };
  // G4double GAr_RIND[20];
  // G4double GAr_ABSL[20];
  // for(int i=0;i<20;i++){
  //   GAr_RIND[i]=GAr_RINDEX;
  //   GAr_ABSL[i]=GAr_ABSORP;
  // }
  //  G4double GAr_ABSL[6]  = { 10000000*m, 10000000*m, 10000000*m, 660*mm, 660*mm, 660*mm}; //atten length
  //  G4double GAr_ABSL[6]  = { 10000000*m, 10000000*m, 10000000*m, 100*mm, 100*mm, 100*mm}; //atten length

  G4double GAr_PP[6]    = { 3.0*eV, 3.0996*eV, 3.1*eV, 9.286*eV, 9.686*eV, 10.186*eV };
  //  G4double GAr_PP[6]    = { 1.9*eV, 2.0*eV, 2.1*eV, 4.0*eV, 6.0*eV, 8.0*eV };
  //  G4double GAr_PP[6]    = { 9.286*eV, 9.686*eV, 10.186*eV };
  G4double GAr_SCINT[6] = {  0.0,  0.000,  0.0,  0.0,  1.0,  0.0 };
  //  G4double GAr_SCINT[6] = {  0.000,  0.000,  0.000,  0.0,  0.000,   0.0 };
  G4double GAr_RIND[6]  = { 1.0004, 1.0004, 1.0004 ,1.0004, 1.0004, 1.0004};
  G4double GAr_ABSL[6]  = { 10000000*m, 10000000*m, 10000000*m, 100000000*m, 10000000*m, 10000000*m}; //atten length
  //  G4double GAr_ABSL[6]  = { 10000000*m, 10000000*m, 10000000*m, 660*mm, 660*mm, 660*mm}; //atten length
  //  G4double GAr_ABSL[6]  = { 10000000*m, 10000000*m, 10000000*m, 100*mm, 100*mm, 100*mm}; //atten length

  //------------------------------------------------------------------------------//

  //--------------------------Liquid Argon Property-------------------------------//

  G4double LAr_PP[103];
  G4double LAr_SCINT[103];
  G4double LAr_RIND[103];
  G4double LAr_RAYL[103];
  G4double LAr_ABSL[103];

  LAr_PP[0] = 3.0995*eV;
  LAr_PP[1] = 3.0996*eV;
  LAr_PP[2] = 3.0997*eV;
  LAr_SCINT[0] = 0.0;
  LAr_SCINT[1] = 1.0;
  LAr_SCINT[2] = 0.0;
  LAr_RIND[0] = 1.233;
  LAr_RIND[1] = 1.233;
  LAr_RIND[2] = 1.233;
  LAr_RAYL[0] = 10000000*m;
  LAr_RAYL[1] = 10000000*m;
  LAr_RAYL[2] = 10000000*m;
  LAr_ABSL[0] = 10000000*m;
  LAr_ABSL[1] = 10000000*m;
  LAr_ABSL[2] = 10000000*m;

  for(G4int i=0;i<100;i++){
    G4double Lambda = (138.-(0.2*(G4double)i))*1e-6;
    LAr_PP[i+3] = LambdaE/Lambda;
    LAr_SCINT[i+3] = GetArScintillationSpectrum(Lambda*1e6);
    LAr_RIND[i+3] = GetLArRefIndex(Lambda);
    LAr_RAYL[i+3] = GetLArRayLength(Lambda,temperature);
    LAr_ABSL[i+3] = GetLArAbsLength(Lambda*1e6);
    //    G4cout << "Lambda = " << Lambda << " LAr_PP = " << LAr_PP[i] << " LAr_SCINTI = " << LAr_SCINT[i] << " LAr_RIND = " << LAr_RIND[i] << " LAr_RAYL = " << LAr_RAYL[i] << " LAr_ABSL = " << LAr_ABSL[i] << G4endl;
  }

  // G4MaterialPropertiesTable *LAr_mt = new G4MaterialPropertiesTable();
  // LAr_mt->AddProperty("FASTCOMPONENT", LAr_PP, LAr_SCINT, 103);
  // LAr_mt->AddProperty("SLOWCOMPONENT", LAr_PP, LAr_SCINT, 103);
  // LAr_mt->AddProperty("RINDEX",        LAr_PP, LAr_RIND,  103);
  // LAr_mt->AddProperty("RAYLEIGH",      LAr_PP, LAr_RAYL,  103);
  // LAr_mt->AddProperty("ABSLENGTH",     LAr_PP, LAr_ABSL,  103);
  // LAr_mt->AddConstProperty("SCINTILLATIONYIELD",0./MeV);
  // LAr_mt->AddConstProperty("RESOLUTIONSCALE",1.0);
  // LAr_mt->AddConstProperty("FASTTIMECONSTANT",7.*ns);
  // LAr_mt->AddConstProperty("SLOWTIMECONSTANT",1600.*ns);
  // LAr_mt->AddConstProperty("YIELDRATIO",0.23);
  // LAr_mat->SetMaterialPropertiesTable(LAr_mt);

  // G4double LAr_PP[6]    = { 3.0*eV, 3.0996*eV, 3.1*eV, 9.286*eV, 9.686*eV, 10.186*eV };
  // G4double LAr_SCINT[6] = { 0.000, 0.000, 0.000, 0.000, 1.000, 0.000 };
  // G4double LAr_RIND[6]  = { 1.233, 1.233, 1.233, 1.356, 1.378, 1.418 };
  // //  G4double LAr_ABSL[6]  = { 10000*m, 10000*m, 10000*m, 10000*m, 10000*m, 10000*m}; // atten length
  // G4double LAr_ABSL[6]  = { 10000*m, 10000*m, 10000*m, 660*mm, 660*mm, 660*mm}; // atten length
  // G4double LAr_RAYL[6]  = {  100*m, 100*m, 100*m, 125.9*cm, 91.07*cm, 56.14*cm}; // Rayleigh scattering length
  // G4MaterialPropertiesTable *LAr_mt = new G4MaterialPropertiesTable();
  // LAr_mt->AddProperty("FASTCOMPONENT", LAr_PP, LAr_SCINT, 6);
  // LAr_mt->AddProperty("SLOWCOMPONENT", LAr_PP, LAr_SCINT, 6);
  // LAr_mt->AddProperty("RINDEX",        LAr_PP, LAr_RIND,  6);
  // LAr_mt->AddProperty("ABSLENGTH",     LAr_PP, LAr_ABSL,  6);
  // LAr_mt->AddProperty("RAYLEIGH",      LAr_PP, LAr_RAYL,  6);
  // LAr_mt->AddConstProperty("SCINTILLATIONYIELD",0./MeV);
  // LAr_mt->AddConstProperty("RESOLUTIONSCALE",1.0);
  // LAr_mt->AddConstProperty("FASTTIMECONSTANT",7.*ns);
  // LAr_mt->AddConstProperty("SLOWTIMECONSTANT",1600.*ns);
  // LAr_mt->AddConstProperty("YIELDRATIO",0.23);
  // LAr_mat->SetMaterialPropertiesTable(LAr_mt);

  //------------------------------------------------------------------------------//

  //------------------------------MEG Fiducial------------------------------------//

  G4MaterialPropertiesTable *MLArFid_mt = new G4MaterialPropertiesTable();
  MLArFid_mt->AddProperty("FASTCOMPONENT", GAr_PP, GAr_SCINT, 6);
  MLArFid_mt->AddProperty("SLOWCOMPONENT", GAr_PP, GAr_SCINT, 6);
  MLArFid_mt->AddProperty("RINDEX",        GAr_PP, GAr_RIND,  6);
  MLArFid_mt->AddProperty("ABSLENGTH",     GAr_PP, GAr_ABSL,  6);
  // MLArFid_mt->AddProperty("FASTCOMPONENT", LAr_PP, LAr_SCINT, 103);
  // MLArFid_mt->AddProperty("SLOWCOMPONENT", LAr_PP, LAr_SCINT, 103);
  // MLArFid_mt->AddProperty("RINDEX",        LAr_PP, LAr_RIND,  103);
  // MLArFid_mt->AddProperty("ABSLENGTH",     LAr_PP, LAr_ABSL,  103);
  // MLArFid_mt->AddProperty("RAYLEIGH",      LAr_PP, LAr_RAYL,  103);
  // MLArFid_mt->AddProperty("FASTCOMPONENT", LAr_PP, LAr_SCINT, 6);
  // MLArFid_mt->AddProperty("SLOWCOMPONENT", LAr_PP, LAr_SCINT, 6);
  // MLArFid_mt->AddProperty("RINDEX",        LAr_PP, LAr_RIND,  6);
  // MLArFid_mt->AddProperty("ABSLENGTH",     LAr_PP, LAr_ABSL,  6);
  // MLArFid_mt->AddProperty("RAYLEIGH",      LAr_PP, LAr_RAYL,  6);
  //   MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",3910./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",3580./MeV); // 1.93 p.e./keVee
  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",500./MeV); // 1.93 p.e./keVee
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",3848.5/MeV); // 5MeV alpha
  // MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",3370./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",300./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",419./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",321./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",572./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",300./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",150.0/MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",196.1/MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",280./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",0./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",978./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",450./MeV);
  //  MLArFid_mt->AddConstProperty("SCINTILLATIONYIELD",302./MeV);
  MLArFid_mt->AddConstProperty("RESOLUTIONSCALE",1.0);
  //  MLArFid_mt->AddConstProperty("RESOLUTIONSCALE",5.);
  MLArFid_mt->AddConstProperty("FASTTIMECONSTANT",7.*ns);
  MLArFid_mt->AddConstProperty("SLOWTIMECONSTANT",1600.*ns);
  MLArFid_mt->AddConstProperty("YIELDRATIO",0.23);
  MLArFid_mat->SetMaterialPropertiesTable(MLArFid_mt);

  //------------------------------------------------------------------------------//

  // Universe - room wall - CONCRETE ****************************************

  //NB: measured INSIDE of lab, therefore have to add twice wall thickness
  G4double wallThick   = 24.*cm;
  G4double worldWidth  = 470.0*cm + 2.*wallThick; // "x"
  G4double worldLength = 690.0*cm + 2.*wallThick; // "y"
  G4double worldHeight = 280.0*cm + 2.*wallThick; // "z"

  G4Box* world_box = new G4Box("world_box", 0.5*worldWidth, 0.5*worldLength, 0.5*worldHeight );
  world_log  = new G4LogicalVolume(world_box, world_mat, "world_log");
  world_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), "world_phys", world_log, NULL, false,0);

  world_log->SetVisAttributes(G4VisAttributes::Invisible);

  G4double labWidth  = worldWidth  - 2.*wallThick; //X
  G4double labLength = worldLength - 2.*wallThick; //Y
  G4double labHeight = worldHeight - 2.*wallThick; //Z

  G4Box* lab_box = new G4Box("lab_box", 0.5*labWidth, 0.5*labLength, 0.5*labHeight );
  lLab  = new G4LogicalVolume(lab_box, lab_mat, "lLab");
  pLab = new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), "pLab", lLab, world_phys, false,0);

  G4VisAttributes* vLab= new G4VisAttributes(white);
  //  lLab->SetVisAttributes(G4VisAttributes::Invisible);
  //  vLab->SetVisibility(true);
  vLab->SetVisibility(false);
  lLab->SetVisAttributes(vLab);

// include room furniture: **************************************************

#include "DMXDetectorRoom.icc"

  //--------Set Visibility-----------
  G4bool IsVisibleCylinder   = false;
  G4bool IsVisibleCollimater = true;
  G4bool IsVisibleHPMT       = true;
  G4bool IsVisibleAcryl      = true;
  G4bool IsVisibleHPMTw      = true;
  G4bool IsVisibleHPMTwls    = true;
  G4bool IsVisibleHGAr       = false;
  G4bool IsVisibleOuterJC    = false;
  G4bool IsVisibleInnerJC    = true;
  G4bool IsVisibleGAr        = true;
  G4bool IsVisibleLAr        = true;
  G4bool IsVisibleFid        = false;
  //---------------------------------

  // Constant Parameter
  G4double PI = 3.14159265359;

  // Detector Geometry

  //****J-Chamber****//
  G4double JCHeight = 1050.0*mm;
  G4double JCThick  = 3.*mm;  
  G4double JCinnerMinR  = 150.*mm;
  G4double JCouterMaxR  = 165.*mm;
  G4double JCinnerMaxR  = JCinnerMinR+JCThick;
  G4double JCouterMinR  = JCouterMaxR-JCThick;
  G4double ChamSPhi = 0.*deg;
  G4double ChamDPhi = 360.*deg;

  //****Gas Argon****//
  G4double GArRadius = 150.0*mm;
  G4double GArHeight = 1050.0*mm;

  //****Liquid Argon****//
  G4double innerLArRadius = 0.*mm;
  G4double outerLArRadius = 150.*mm;
  G4double HeightLAr = 1050.*mm;
  G4double startAngleLAr  = 0.*deg;
  G4double spanningAngleLAr = 360.*deg;

  //****PMT****//
  G4double MPMTRadius = 27.8*mm; // MEG
  //  G4double MPMTRadius = 36.*(1.+(2./sqrt(3.)))*mm; // Triple J
  //  G4double MPMTRadius = 25.4*mm; // H1161
  //  G4double MPMTActiveRadius = 23.0*mm; // H1161
  G4double MPMTActiveRadius = 27.8*mm; // MEG
  //  G4double MPMTActiveRadius = 25.4*mm; // MEG
  //  G4double MPMTActiveRadius = 27.8*mm; // MEG
  //  G4double MPMTActiveRadius = 36.*(1.+(2./sqrt(3.)))*mm; // TripleJ 
  G4double MPMTHeight = 1.0*mm;

  //****PMT window****//
  G4double MPMTwThick = 0.1*mm; // MEG
  G4double MPMTwRadius = MPMTRadius; // MEG

  //****PMT window WLS****//
  G4double MPwlsThick = 0.001*mm; // MEG
  //  G4double MPwlsThick = 0.*mm; // MEG
  G4double MPwlsRadius = MPMTRadius; // MEG

  //****fuducial****//
  G4double MFidRadius = MPMTRadius;
  G4double MFidHeight = 60.0*mm; // RunV
  //  G4double MFidHeight = 55.0*mm; // RunV
  //  G4double MFidHeight = 80.0*mm; // RunV
  //  G4double MFidHeight = 100.0*mm; // RunV
  //  G4double MFidHeight = 60.0*4*mm; // RunV
  //  G4double MFidHeight = 80.0*mm; // Gas Setup

  //***************** Reflector WLS Section ************************//
  //----------------- MEG Refector WLS----------------//
  G4double MrefwlsThick = 0.001*mm;
  G4double MrefwlsRadius = MFidRadius;
  G4double MrefwlsHeight = MFidHeight;

  //***************** Reflector Section *************************//
  //----------------- MEG Refector -------------------//
  G4double MReflectorThick = 0.1*mm;
  //  G4double MReflectorRadius = MFidRadius+MrefwlsThick;
  G4double MReflectorRadius = MFidRadius;
  //  G4double MReflectorHeight = MFidHeight+MrefwlsThick*2;
  G4double MReflectorHeight = MFidHeight;

  //***************** Outer Reflector Section ******************//

  //******************MEG Cylinder***********************//
  //  G4double MPMTCylinder_Thick = 0.1*mm;
  G4double MPMTCylinder_Thick = (MPMTRadius-MPMTActiveRadius)+MrefwlsThick+MReflectorThick;
  G4double MPMTCylinder_TubeHeight = MPMTHeight;
  G4double MPMTCylinder_PlateHeight = 32.4*mm;
  //  G4double MPMTCylinder_PlateHeight = 0.1*mm;

  //*****************VETO PMT FIX PLATE********************//

  G4double vetofixplatethick=3.0*mm;
  G4double vetofixplateradius=133.*mm;
  G4double vetofixplateVar=189./2*mm;

  //******************MEG Detector***********************//
  G4double MEGXDis = 30.*mm; // RunV
  //  G4double MEGXDis = 40.*mm; // RunV
  //  G4double MEGXDis = 27.5*mm; // RunV
  //  G4double MEGXDis = 30.*4*mm; // RunV
  //  G4double MEGXDis = 40.*mm; // Gas Setup
  //  G4double MEGYDis = 84.12*mm;
  //  G4double MEGZDis = vetofixplateVar-2.5-MPMTRadius;  // 64.2
  G4double MEGYDis = 0.*mm;
  G4double MEGZDis = 0.*mm;  // 64.2

  //******************VETO MEG***********************//

  G4double vetoRadius = 27.8*mm;
  G4double vetowlsThick = 0.001*mm;
  G4double vetowinThick = 1.0*mm;
  G4double vetopmtThick = 1.0*mm;
  G4double vetoDis = 95.*mm;
  G4double vetoVar = vetofixplateVar-MPwlsThick-MPMTwThick-MPMTCylinder_TubeHeight-MPMTCylinder_PlateHeight;

  // Rotation

  G4double rotA = 90.;
  G4RotationMatrix* rotP = new G4RotationMatrix();
  rotP->rotateY(rotA*deg); 

  G4double rotB = 90.;
  G4RotationMatrix* rotQ = new G4RotationMatrix();
  rotQ->rotateY(-rotB*deg); 

  G4RotationMatrix* rotJ = new G4RotationMatrix();
  rotJ->rotateY(180.*deg); 

  // Now start with detector assembly:

  //*******************************************************//

  //*******************J-Chamber**************************// 

  G4Tubs* JCinnerTube = new G4Tubs("JCinnerTube",JCinnerMinR,JCinnerMaxR,0.5*JCHeight,0.*deg,360.*deg);
  G4Tubs* JCouterTube = new G4Tubs("JCouterTube",JCouterMinR,JCouterMaxR,0.5*JCHeight,0.*deg,360.*deg);

  G4UnionSolid* JCTube = new G4UnionSolid("JCTube",JCinnerTube,JCouterTube,0,G4ThreeVector(0.,0.,0.));

  lJC  = new G4LogicalVolume(JCTube, JC_mat, "lJC");
  pJC = new G4PVPlacement(0,G4ThreeVector(0.,0.,0.),lJC,"pJC", lLab, false,0);

  G4VisAttributes* vJC = new G4VisAttributes(white);
  //  vJC->SetVisibility(true);
  vJC->SetVisibility(false);
  lJC->SetVisAttributes(vJC);

  //*******************************************************//

  //*******************GAr Tube***************************//

  G4Tubs* GAr_tube = new G4Tubs("GAr_tube",0,GArRadius,0.5*GArHeight, 0.*deg, 360.*deg);

  lHGAr = new G4LogicalVolume(GAr_tube, GAr_mat, "lHGAr");
  //  pHGAr = new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), lHGAr,"pHGAr",lLab, false,0);

  G4VisAttributes* vHGAr = new G4VisAttributes(white);
  //  vHGAr->SetVisibility(IsVisibleHGAr);
  vHGAr->SetVisibility(false);
  lHGAr->SetVisAttributes(vHGAr);
  //******************************************************//

  //*******************LAr Tube***************************//

  G4Tubs* LAr_tube=new G4Tubs("LAr_tube",innerLArRadius,outerLArRadius,0.5*HeightLAr,startAngleLAr,spanningAngleLAr);
  lLAr=new G4LogicalVolume(LAr_tube,LAr_mat,"lLAr"); 
  pLAr=new G4PVPlacement(0,G4ThreeVector(0.,0.,0.),lLAr,"pLAr",lLab,false,0);

  G4VisAttributes* vLAr=new G4VisAttributes(blue);
  vLAr->SetVisibility(false);
  lLAr->SetVisAttributes(vLAr);
  //******************************************************//

  //**************MEG Reflector Substruct Tube****************//

  G4double MRSubR = 5.*mm;
  G4double MRSubH = 10.*mm;

  G4Tubs* MRSub_tube=new G4Tubs("MRSub_tube",0,MRSubR,0.5*MRSubH,0.*deg, 360.*deg);

  //**************MEG Reflector Substruct Box****************//

  G4double MRSubXYZ = 1.*mm;
  
  G4Box* MRSub_Box=new G4Box("MRSub_Box",MRSubXYZ*0.5,MRSubXYZ*0.5,MRSubXYZ*0.5);

  //**************MEG Reflector Section****************//

  G4Tubs* MReflector_tube = new G4Tubs("MReflector_tube",MReflectorRadius,MReflectorRadius+MReflectorThick,
				       0.5*MReflectorHeight, 0.*deg, 360.*deg);

  G4SubtractionSolid* MReflector_tube_sub = new G4SubtractionSolid("MReflector_tube_sub",
   								   MReflector_tube,MRSub_tube,rotP,
   								   G4ThreeVector(MReflectorRadius,0,0));

  // G4SubtractionSolid* MReflector_tube_sub = new G4SubtractionSolid("MReflector_tube_sub",
  // 								   MReflector_tube,MRSub_Box,rotP,
  // 								   G4ThreeVector(-MReflectorRadius,0,0));
  
  lMReflector = new G4LogicalVolume(MReflector_tube, Ref_mat, "lMReflector");
  //  lMReflector = new G4LogicalVolume(MReflector_tube_sub, Ref_mat, "lMReflector");
  pMReflector = new G4PVPlacement(rotP, G4ThreeVector(0.,-MEGYDis,-MEGZDis),lMReflector,"pMReflector",lLAr, false,0);

  G4VisAttributes* vMReflector=new G4VisAttributes(white);
  vMReflector->SetVisibility(true);
  lMReflector->SetVisAttributes(vMReflector);

  //******************************************************//

  //***************MEG Reflector WLS Secssion*******************//

  G4Tubs* Mrefwls_tube = new G4Tubs("Mrefwls_tube",MrefwlsRadius,MrefwlsRadius+MrefwlsThick,
				     0.5*MrefwlsHeight, 0.*deg, 360.*deg);  

  G4SubtractionSolid* Mrefwls_tube_sub = new G4SubtractionSolid("Mrefwls_tube_sub",
  								Mrefwls_tube,MRSub_tube,rotP,
  								G4ThreeVector(MrefwlsRadius,0,0));

  // G4SubtractionSolid* Mrefwls_tube_sub = new G4SubtractionSolid("Mrefwls_tube_sub",
  // 								Mrefwls_tube,MRSub_Box,rotP,
  // 								G4ThreeVector(-MrefwlsRadius,0,0));

  //  lMrefwls = new G4LogicalVolume(Mrefwls_tube, wls_mat, "lMrefwls");
  lMrefwls = new G4LogicalVolume(Mrefwls_tube, RefWLS_mat, "lMrefwls");
  //  lMrefwls = new G4LogicalVolume(Mrefwls_tube_sub, wls_mat, "lMrefwls");
  //  pMrefwls = new G4PVPlacement(rotP, G4ThreeVector(0.,-MEGYDis,-MEGZDis),lMrefwls,"pMrefwls",lLAr, false,0);

  G4VisAttributes* vMrefwls=new G4VisAttributes(blue);
  vMrefwls->SetVisibility(false);
  lMrefwls->SetVisAttributes(vMrefwls);

  //******************************************************//

  //*****************MEG PMT SECSSION*********************// 

  //----------------MEGPMT WLS---------------------//
  G4Tubs* Mwls_plate = new G4Tubs("Mwls_plate",0.*cm, MPwlsRadius, 0.5*MPwlsThick, 0.*deg, 360.*deg);
  lMPwls1 = new G4LogicalVolume(Mwls_plate, wls_mat, "lMPwls1");
  lMPwls2 = new G4LogicalVolume(Mwls_plate, wls_mat, "lMPwls2");

  //--------------MEGPMT Window--------------------//
  G4Tubs* MPMTw_plate = new G4Tubs("MPMTw1_plate",0.*cm, MPMTwRadius, 0.5*MPMTwThick, 0.*deg, 360.*deg);
  lMPMTw = new G4LogicalVolume(MPMTw_plate, PMTS_mat, "lMPMTw");

  //---------------MEG PMT-------------------------//

  //  G4Tubs* MPMT_plate = new G4Tubs("MPMT_plate",0.*cm, MPMTRadius, 0.5*MPMTHeight, 0.*deg, 360.*deg);
  G4Tubs* MPMT_plate = new G4Tubs("MPMT_plate",0.*cm, MPMTActiveRadius, 0.5*MPMTHeight, 0.*deg, 360.*deg);
  lMPMT = new G4LogicalVolume(MPMT_plate, SP_mat, "lMPMT");

  //---------------MEG PMT Cylinder----------------//

  //  G4Tubs* MPMTCylinder_tube = new G4Tubs("MPMTCylinder_tube",MPMTRadius, MPMTRadius+MPMTCylinder_Thick, 0.5*MPMTCylinder_TubeHeight, 0.*deg, 360.*deg);
  G4Tubs* MPMTCylinder_tube = new G4Tubs("MPMTCylinder_tube",MPMTActiveRadius, MPMTActiveRadius+MPMTCylinder_Thick, 0.5*MPMTCylinder_TubeHeight, 0.*deg, 360.*deg);
  //  G4Tubs* MPMTCylinder_plate = new G4Tubs("MPMTCylinder_plate", 0, MPMTRadius+MPMTCylinder_Thick, 0.5*MPMTCylinder_PlateHeight, 0.*deg, 360.*deg);
  G4Tubs* MPMTCylinder_plate = new G4Tubs("MPMTCylinder_plate", 0, MPMTActiveRadius+MPMTCylinder_Thick, 0.5*MPMTCylinder_PlateHeight, 0.*deg, 360.*deg);
  G4UnionSolid* MPMTCylinder = new G4UnionSolid("MPMTCylinder",MPMTCylinder_tube,MPMTCylinder_plate,0,G4ThreeVector(0,0,0.5*MPMTCylinder_TubeHeight+0.5*MPMTCylinder_PlateHeight));
  lMPMTCylinder = new G4LogicalVolume(MPMTCylinder, SP_mat, "lMPMTCylinder");

  //****************MEG Detector*********************// 

  //*************MPMT Fiducial Tube*****************//

  G4Tubs* MFid_tube=new G4Tubs("MFid_tube",0,MFidRadius,0.5*MFidHeight,0.*deg,360.*deg);
  lMFid=new G4LogicalVolume(MFid_tube,MLArFid_mat,"lMFid"); 
  pMFid=new G4PVPlacement(rotP,G4ThreeVector(0.,-MEGYDis,-MEGZDis),lMFid,"pMFid",lLAr,false,0);

  G4VisAttributes* vMFid=new G4VisAttributes(blue);
  vMFid->SetVisibility(false);
  lMFid->SetVisAttributes(vMFid);

  //************************************************//

  pMPwls1 = new G4PVPlacement(rotQ,G4ThreeVector(MEGXDis+MPwlsThick*0.5,-MEGYDis,-MEGZDis),
      			      lMPwls1, "pMPwls1", lLAr, true, 0);
  pMPwls2 = new G4PVPlacement(rotP,G4ThreeVector(-MEGXDis-MPwlsThick*0.5,-MEGYDis,-MEGZDis),
  			      lMPwls2, "pMPwls2", lLAr, true, 0);
  pMPMTw1 = new G4PVPlacement(rotQ,G4ThreeVector(MEGXDis+MPwlsThick+MPMTwThick*0.5,-MEGYDis,-MEGZDis),
      			      lMPMTw, "pMPMTw1", lLAr, true, 0);
  //  pMPMTw1 = new G4PVPlacement(rotQ,G4ThreeVector(MEGXDis+MPMTwThick*0.5,-MEGYDis,-MEGZDis),
  //    			      lMPMTw, "pMPMTw1", lLAr, true, 0);
  pMPMTw2 = new G4PVPlacement(rotP,G4ThreeVector(-MEGXDis-MPwlsThick-MPMTwThick*0.5,-MEGYDis,-MEGZDis),
			      lMPMTw, "pMPMTw2", lLAr, true, 0); 
  //  pMPMTw2 = new G4PVPlacement(rotP,G4ThreeVector(-MEGXDis-MPMTwThick*0.5,-MEGYDis,-MEGZDis),
  //      			      lMPMTw, "pMPMTw2", lLAr, true, 0); 
  pMPMT1 = new G4PVPlacement(rotQ,G4ThreeVector(MEGXDis+MPwlsThick+MPMTwThick+MPMTHeight*0.5,-MEGYDis,-MEGZDis),
  			     lMPMT, "pMPMT1", lLAr, true, 0); 
  //  pMPMT1 = new G4PVPlacement(rotQ,G4ThreeVector(MEGXDis+MPMTwThick+MPMTHeight*0.5,-MEGYDis,-MEGZDis),
  //			     lMPMT, "pMPMT1", lLAr, true, 0); 
  pMPMT2 = new G4PVPlacement(rotP,G4ThreeVector(-MEGXDis-MPwlsThick-MPMTwThick-MPMTHeight*0.5,-MEGYDis,-MEGZDis),
        			     lMPMT, "pMPMT2", lLAr, true, 0); 
  //  pMPMT2 = new G4PVPlacement(rotP,G4ThreeVector(-MEGXDis-MPMTwThick-MPMTHeight*0.5,-MEGYDis,-MEGZDis),
  //      			     lMPMT, "pMPMT2", lLAr, true, 0); 
  pMPMTCylinder1 = new G4PVPlacement(rotQ,G4ThreeVector(MEGXDis+MPwlsThick+MPMTwThick+MPMTCylinder_TubeHeight*0.5,-MEGYDis,-MEGZDis),
    				     lMPMTCylinder, "pMPMTCylinder1", lLAr, true, 0); 
  //  pMPMTCylinder1 = new G4PVPlacement(rotQ,G4ThreeVector(MEGXDis+MPMTwThick+MPMTCylinder_TubeHeight*0.5,-MEGYDis,-MEGZDis),
  //  				     lMPMTCylinder, "pMPMTCylinder1", lLAr, true, 0); 
  pMPMTCylinder2 = new G4PVPlacement(rotP,G4ThreeVector(-MEGXDis-MPwlsThick-MPMTwThick-MPMTCylinder_TubeHeight*0.5,-MEGYDis,-MEGZDis),
				     lMPMTCylinder, "pMPMTCylinder2", lLAr, true, 0); 
  //  pMPMTCylinder2 = new G4PVPlacement(rotP,G4ThreeVector(-MEGXDis-MPMTwThick-MPMTCylinder_TubeHeight*0.5,-MEGYDis,-MEGZDis),
  //      				     lMPMTCylinder, "pMPMTCylinder2", lLAr, true, 0); 

  G4VisAttributes* vMPwls = new G4VisAttributes(purple);
  vMPwls->SetForceSolid(true);
  vMPwls->SetVisibility(true); 
  lMPwls1->SetVisAttributes(vMPwls);  
  lMPwls2->SetVisAttributes(vMPwls);  

  G4VisAttributes* vMPMTw = new G4VisAttributes(yellow);
  vMPMTw->SetForceSolid(true);
  vMPMTw->SetVisibility(true);
  lMPMTw->SetVisAttributes(vMPMTw);  

  G4VisAttributes* vMPMT = new G4VisAttributes(cyan);
  vMPMT->SetForceSolid(true);
  vMPMT->SetVisibility(true);
  lMPMT->SetVisAttributes(vMPMT);  

  G4VisAttributes* vMPMTCylinder = new G4VisAttributes(white);
  vMPMTCylinder->SetForceSolid(true);
  vMPMTCylinder->SetVisibility(true); 
  lMPMTCylinder->SetVisAttributes(vMPMTCylinder);  

  //*****************VETO PMT FIX PLATE********************//

  G4Tubs* vetofixplate = new G4Tubs("vetoCylinder_tube",0., vetofixplateradius, 0.5*vetofixplatethick, 0.*deg, 360.*deg);

  lvetofixplate = new G4LogicalVolume(vetofixplate, SP_mat, "lvetofixplate");
  //  pvetofixplate1 = new G4PVPlacement(0, G4ThreeVector(0.,0.,vetofixplateVar+0.5*vetofixplatethick),
  //  				     lvetofixplate, "pvetofixplate1", lLAr, false,0);
  //  pvetofixplate2 = new G4PVPlacement(0, G4ThreeVector(0.,0.,-vetofixplateVar-0.5*vetofixplatethick),
  //  				     lvetofixplate, "pvetofixplate2", lLAr, false,0);

  G4VisAttributes* vvetofixplate = new G4VisAttributes(white);
  vvetofixplate->SetForceSolid(false);
  vvetofixplate->SetVisibility(false); 
  lvetofixplate->SetVisAttributes(vvetofixplate);  

  //******************************************************//

  //********Optical Condition*********//

  //***********PMT WLS && Window****************//

  G4OpticalSurface * OpMPMTWSurface = new G4OpticalSurface
    ("OpMPMTWSurface", unified, polished, dielectric_dielectric, 0.01*deg);

  G4LogicalSkinSurface *MPMTwSurface = new G4LogicalSkinSurface("MPMTwSurface",lMPMTw,OpMPMTWSurface);

  G4OpticalSurface * OpMWLSurface2 = new G4OpticalSurface
    //    ("OpMWLSurface", unified, polished, dielectric_dielectric, 0.01*deg);
    ("OpWLSurface2", unified, ground, dielectric_dielectric, 89.99*deg);
    //    ("OpWLSurface", unified, ground, dielectric_dielectric, 0.01*deg);

  G4OpticalSurface * OpMWLSurface1 = new G4OpticalSurface
    ("OpMWLSurface1", unified, polished, dielectric_dielectric, 0.01*deg);
    //    ("OpWLSurface", unified, ground, dielectric_dielectric, 89.99*deg);
    //    ("OpWLSurface", unified, ground, dielectric_dielectric, 0.01*deg);

  G4MaterialPropertiesTable *MPMTwls_mt = new G4MaterialPropertiesTable();
  G4double MPMTwls_PP[4]   = { 2.0*eV, 3.5*eV, 9.0*eV, 11.0*eV };
  //  G4double MPMTwls_RIND[4] = {   1.486,   1.486,   1.486,    1.486 };     // ref index
  G4double MPMTwls_RIND[4] = {   1.635,   1.635,   1.635,    1.635 };     // ref index
  //  G4double MPMTwls_RIND[4] = {   1.233,   1.233,   1.233,    1.233 };     // ref index
  G4double MPMTwls_SPIKE[4] = {  0.0,   0.0,    0.0,     0.0 };     // ref index
  G4double MPMTwls_LOBE[4] = {  1.0,   1.0,    0.0,     0.0 };     // ref index
  MPMTwls_mt->AddProperty("RINDEX", MPMTwls_PP, MPMTwls_RIND, 4);
  MPMTwls_mt->AddProperty("SPECULARSPIKECONSTANT", MPMTwls_PP, MPMTwls_SPIKE, 4);
  MPMTwls_mt->AddProperty("SPECULARLOBECONSTANT", MPMTwls_PP, MPMTwls_LOBE, 4);
  OpMWLSurface2->SetMaterialPropertiesTable(MPMTwls_mt);

  const G4int Wls_NUMENTRIES = 3;
  G4double Wls_PP[Wls_NUMENTRIES]   = { 2.0*eV, 3.5*eV ,11.0*eV};
  G4double Wls_RIND[Wls_NUMENTRIES] = { 1.486, 1.486 ,1.486};     // ref index
  //  G4double Wls_RIND[Wls_NUMENTRIES] = { 1.635, 1.635 ,1.635};     // ref index
  //  G4double Wls_RIND[Wls_NUMENTRIES] = { 1.233, 1.233 ,1.233};     // ref index
  //  G4double Wls_RIND[Wls_NUMENTRIES] = { 1.0004, 1.0004 ,1.0004};     // ref index
  G4MaterialPropertiesTable *Wls_mt = new G4MaterialPropertiesTable();
  Wls_mt->AddProperty("RINDEX", Wls_PP, Wls_RIND, Wls_NUMENTRIES);
  OpMWLSurface1->SetMaterialPropertiesTable(Wls_mt);

  G4LogicalSkinSurface *MPwlsSurface1 = new G4LogicalSkinSurface("MPwlsSurface1",lMPwls1,OpMWLSurface2);
  G4LogicalSkinSurface *MPwlsSurface2 = new G4LogicalSkinSurface("MPwlsSurface2",lMPwls2,OpMWLSurface2);

  //******************************************************//

  //*************Reflector WLS && Reflector***************//

  //------------------ MRefeflector Surface -------------------------
  G4OpticalSurface * OpMRefSurface = new G4OpticalSurface
    ("OpMRefSurface", unified, polished, dielectric_metal, 0.01*deg);
    //    ("OpMRefSurface", unified, ground, dielectric_metal, 90.0*deg);

  G4double MRef_PP[4]   = { 2.0*eV, 3.5*eV, 9.0*eV, 11.0*eV };
  G4double MRef_REFL[4] = {   0.00,   0.00,    0.0,     0.0 };     // ref index
  //  G4double MRef_REFL[4] = {   0.99,   0.99,    0.0,     0.0 };     // ref index
  //  G4double MRef_REFL[4] = {   0.98,   0.98,    0.0,     0.0 };     // ref index
  //  G4double MRef_REFL[4] = {   1.00,   1.00,    0.0,     0.0 };     // ref index
  //  G4double MRef_REFL[4] = {   0.99,   0.99,    0.0,     0.0 };     // ref index
  //  G4double MRef_REFL[4] = {   0.9,   0.9,    0.0,     0.0 };     // ref index
  //  G4double MRef_REFL[4] = {   0.95,   0.95,    0.0,     0.0 };     // ref index
  //  G4double MRef_SPIKE[4] = {0.9, 0.9, 0.0, 0.0};
  //  G4double MRef_LOBE[4] = {0.1, 0.1, 0.0, 0.0};
  G4double MRef_SPIKE[4] = {1.0, 1.0, 0.0, 0.0};
  G4double MRef_LOBE[4] = {0.0, 0.0, 1.0, 1.0};

  G4MaterialPropertiesTable *MRef_mt = new G4MaterialPropertiesTable();
  MRef_mt->AddProperty("REFLECTIVITY", MRef_PP, MRef_REFL, 4);
  //  MRef_mt->AddProperty("SPECULARSPIKECONSTANT",MRef_PP,MRef_SPIKE, 4);
  //  MRef_mt->AddProperty("SPECULARLOBECONSTANT",MRef_PP,MRef_LOBE, 4);
  OpMRefSurface->SetMaterialPropertiesTable(MRef_mt);

  G4LogicalBorderSurface *MRefSurface;
  //  MRefSurface = new G4LogicalBorderSurface("MRefSurface1_1", pMrefwls, pMReflector, OpMRefSurface);
  MRefSurface = new G4LogicalBorderSurface("MRefSurface1_1", pMFid, pMReflector, OpMRefSurface);

  //------------------ MReflector WLS Surface -------------------------
  G4OpticalSurface *OpMRefwlsSurface = new G4OpticalSurface
    //    ("OpMRefwlsSurface", unified, polished, dielectric_dielectric, 0.01*deg);
    //    ("OpMRefwlsSurface", unified, ground, dielectric_dielectric, 10.5*deg);
    ("OpMRefwlsSurface", unified, ground, dielectric_dielectric, 89.99*deg);
    //    ("OpMRefwlsSurface", unified, ground, dielectric_dielectric, 2.0*deg);
    //    ("OpMRefwlsSurface", unified, ground, dielectric_dielectric, 5.0*deg);

  G4MaterialPropertiesTable *MRefwls_mt = new G4MaterialPropertiesTable();
  G4double MRefwls_PP[4]   = { 2.0*eV, 3.5*eV, 9.0*eV, 11.0*eV };
  //  G4double MRefwls_RIND[4] = {   1.486,   1.486,   1.486,    1.486 };     // ref index
  //  G4double MRefwls_RIND[4] = {   1.0004,   1.0004,   1.0004,    1.0004 };     // ref index
  G4double MRefwls_RIND[4] = {   1.635,   1.635,   1.635,    1.635 };     // ref index
  //  G4double MRefwls_RIND[4] = {   1.233,   1.233,   1.233,    1.233 };     // ref index
  G4double MRefwls_SPIKE[4] = {  0.0,   0.0,    0.0,     0.0 };     // ref index
  G4double MRefwls_LOBE[4] = {  1.0,   1.0,    0.0,     0.0 };     // ref index
  MRefwls_mt->AddProperty("RINDEX", MRefwls_PP, MRefwls_RIND, 4);
  MRefwls_mt->AddProperty("SPECULARSPIKECONSTANT", MRefwls_PP, MRefwls_SPIKE, 4);
  MRefwls_mt->AddProperty("SPECULARLOBECONSTANT", MRefwls_PP, MRefwls_LOBE, 4);
  OpMRefwlsSurface->SetMaterialPropertiesTable(MRefwls_mt);

  G4LogicalSkinSurface *MRefwlsSurface = new G4LogicalSkinSurface("MRefwlsSurface",lMrefwls,OpMRefwlsSurface);

  //******************************************************//

  //***************STOP TOTAL REFLECTION******************//

  G4double strinnerRadius = MFidRadius;
  G4double strouterRadius = MFidRadius + MrefwlsThick;
  G4double strHeight1 = MPwlsThick+MPMTwThick;

  G4Tubs* str_tube1 = new G4Tubs("str_tube1", strinnerRadius, strouterRadius, 0.5*strHeight1, 0.*deg, 360.*deg);
  G4Tubs* str_tube2 = new G4Tubs("str_tube2", strinnerRadius, strouterRadius, 0.5*strHeight1, 0.*deg, 360.*deg);

  lstr1 = new G4LogicalVolume(str_tube1, str_mat, "lstr1");
  lstr2 = new G4LogicalVolume(str_tube2, str_mat, "lstr2");
  pstr1 = new G4PVPlacement(rotP, G4ThreeVector(-(MEGXDis+(MPwlsThick+MPMTwThick)*0.5),-MEGYDis,-MEGZDis),
        			    lstr1, "pstr1", lLAr, false, 0);
  //  pstr1 = new G4PVPlacement(rotP, G4ThreeVector(-(MEGXDis+(MPMTwThick)*0.5),-MEGYDis,-MEGZDis),
  //      			    lstr1, "pstr1", lLAr, false, 0);
  //  pstr11 = new G4PVPlacement(0, G4ThreeVector(0,0,distance+PwlsThick+JPMTwHeight*0.5),
  //			     lstr2, "pstr11", lLAr, false, 0);
  pstr2 = new G4PVPlacement(rotP, G4ThreeVector(MEGXDis+(MPwlsThick+MPMTwThick)*0.5,-MEGYDis,-MEGZDis),
    			    lstr2, "pstr2", lLAr, false, 0);
  //  pstr2 = new G4PVPlacement(rotP, G4ThreeVector(MEGXDis+(MPMTwThick)*0.5,-MEGYDis,-MEGZDis),
  //  			    lstr2, "pstr2", lLAr, false, 0);
  //  pstr22 = new G4PVPlacement(0, G4ThreeVector(0,0,-distance-PwlsThick-JPMTwHeight*0.5),
  //  			     lstr2, "pstr22", lLAr, false, 0);

  G4VisAttributes* vstr = new G4VisAttributes(white);
  vstr->SetVisibility(true);
  vstr->SetForceSolid(true);
  lstr1->SetVisAttributes(vstr);
  lstr2->SetVisAttributes(vstr);

  //******************************************************//

  //******************************************************//

  // Rotation

  G4double rotAA = 90.;
  G4RotationMatrix* rotPP = new G4RotationMatrix();
  rotPP->rotateY(-rotAA*deg); 

  //*****************Capton Sheet*************************//

  G4double capThick        = 0.06*mm;
  G4double capXY           = 2.0*mm;
  //  G4double capdistance     = MPMTRadius-3.2*mm;
  //  G4double capdistance     = MPMTRadius-pulThick-0.5*capThick;
  G4double capdistance     = MPMTRadius-0.1-0.5*capThick;

  G4Box* Cap_Box=new G4Box("Cap_Box",capXY*0.5,capXY*0.5,capThick*0.5);

  lcap = new G4LogicalVolume(Cap_Box, Cap_mat, "lcap");
  //  pcap = new G4PVPlacement(rotPP,G4ThreeVector(-capdistance,0.,0.), lcap, "pcap", lMFid, false, 0);

  G4VisAttributes* vcap = new G4VisAttributes(cyan);
  vcap->SetVisibility(true);
  vcap->SetForceSolid(true);
  lcap->SetVisAttributes(vcap);

  //******************************************************//  

  G4double pulThick        = 0.1*mm;
  //  G4double colHeight       = 5.0*mm;
  G4double colHeight       = 8.0*mm;

  //*****************ESR Sheet*************************//

  //  G4double ESRThick        = 0.018*mm;
  G4double ESRThick        = 0.015*mm; // nantokaSOU
  G4double ESRXY           = 2.0*mm;
  //  G4double ESRdistance     = MPMTRadius-3.2*mm;
  G4double ESRdistance     = MPMTRadius-pulThick-0.5*ESRThick;
  //  G4double ESRdistance     = MPMTRadius-0.1-0.5*ESRThick;

  G4Box* ESR_Box=new G4Box("ESR_Box",ESRXY*0.5,ESRXY*0.5,ESRThick*0.5);
  //  G4Box* ESR_Box=new G4Box("ESR_Box",ESRXY*0.25,ESRXY*0.5,ESRThick*0.5);

  lESR = new G4LogicalVolume(ESR_Box, Ref_mat, "lESR");
  //  pESR = new G4PVPlacement(rotPP,G4ThreeVector(-ESRdistance,0.,0.), lESR, "pESR", lMFid, false, 0);
  //  pESR = new G4PVPlacement(rotPP,G4ThreeVector(-ESRdistance,0.,0.5), lESR, "pESR", lMFid, false, 0);

  G4VisAttributes* vESR = new G4VisAttributes(cyan);
  vESR->SetVisibility(true);
  vESR->SetForceSolid(true);
  lESR->SetVisAttributes(vESR);

  //****************************************************//  

  //*****************Lap Sheet*************************//

  G4double LapThick        = 0.011*mm;
  G4double LapXY           = 4.0*mm;
  //  G4double Lapdistance     = MPMTRadius-pulThick-colHeight-0.5*LapThick;
  G4double Lapdistance     = MPMTRadius-pulThick-0.5*LapThick;

  G4Box* Lap_Box=new G4Box("Lap_Box",LapXY*0.5,LapXY*0.5,LapThick*0.5);

  lLap = new G4LogicalVolume(Lap_Box, Ref_mat, "lLap");
  //  pLap = new G4PVPlacement(rotPP,G4ThreeVector(-Lapdistance,0.,0.), lLap, "pLap", lMFid, false, 0);

  G4VisAttributes* vLap = new G4VisAttributes(cyan);
  vLap->SetVisibility(true);
  vLap->SetForceSolid(true);
  lLap->SetVisAttributes(vLap);

  //***************************************************//  

  //*****************Paraloid Sheet*************************//

  G4double ParaloidThick        = 0.036*mm;
  G4double ParaloidXY           = 4.0*mm;
  G4double Paraloiddistance     = MPMTRadius-pulThick-0.5*ParaloidThick;

  G4Box* Paraloid_Box=new G4Box("Paraloid_Box",ParaloidXY*0.5,ParaloidXY*0.5,ParaloidThick*0.5);

  lParaloid = new G4LogicalVolume(Paraloid_Box, Paraloid_mat, "lParaloid");
  //  pParaloid = new G4PVPlacement(rotPP,G4ThreeVector(-Paraloiddistance,0.,0.), lParaloid, "pParaloid", lMFid, false, 0);

  G4VisAttributes* vParaloid = new G4VisAttributes(cyan);
  vParaloid->SetVisibility(true);
  vParaloid->SetForceSolid(true);
  lParaloid->SetVisAttributes(vParaloid);

  //***************************************************//  

  //*****************Am241 Sheet*************************//

  G4double Am241Thick        = 0.1*mm;
  G4double Am241XY           = 4.0*mm;
  G4double Am241distance     = MPMTRadius-pulThick-0.5*Am241Thick;

  G4Box* Am241_Box=new G4Box("Am241_Box",Am241XY*0.5,Am241XY*0.5,Am241Thick*0.5);

  lAm241 = new G4LogicalVolume(Am241_Box, americium_mat, "lAm241");
  //  pAm241 = new G4PVPlacement(rotPP,G4ThreeVector(-Am241distance,0.,0.), lAm241, "pAm241", lMFid, false, 0);

  G4VisAttributes* vAm241 = new G4VisAttributes(cyan);
  vAm241->SetVisibility(true);
  vAm241->SetForceSolid(true);
  lAm241->SetVisAttributes(vAm241);

  //***************************************************//  
  
  //******************Collimater**************************//
 
  G4double colThick    = 1.0*mm; // with collimeter for alpha without ref
  G4double colpThick   = 0.1*mm; // with collimeter for alpha without ref
  G4double colRadius   = 3.0*mm; // with collimeter for alpha without ref
  //  G4double colHeight   = 8.0*mm; // with collimeter for alpha without ref
  //  G4double coldistance = MPMTRadius-pulThick-0.5*colHeight;
  //  G4double coldistance = MPMTRadius-pulThick-ESRThick-0.5*colHeight;
  //  G4double coldistance = MPMTRadius-colpThick-LapThick-0.5*colHeight;
  G4double coldistance = MPMTRadius-colpThick-0.5*colHeight;

  G4Tubs* col_tube = new G4Tubs("col_tube", colRadius-colThick*0.5, colRadius+colThick*0.5, 
                                0.5*colHeight, 0.*deg, 360.*deg); 
  G4Tubs* col_plate = new G4Tubs("col_plate", 0., colRadius+colThick*0.5, 
                                0.5*colpThick, 0.*deg, 360.*deg); 

  G4UnionSolid* col = new G4UnionSolid("col",col_tube,col_plate,0,
				       G4ThreeVector(0,0,-0.5*colHeight-0.5*colpThick));

  lcol = new G4LogicalVolume(col, col_mat, "lcol");
  //  lcol = new G4LogicalVolume(col_tube, col_mat, "lcol");
  pcol = new G4PVPlacement(rotPP,G4ThreeVector(-coldistance,0.,0.), lcol, "pcol", lMFid, false, 0);
  //  pcol = new G4PVPlacement(rotPP,G4ThreeVector(-coldistance,0.,5.0), lcol, "pcol", lMFid, false, 0);
  
  G4VisAttributes* vcol = new G4VisAttributes(red);
  vcol->SetVisibility(true);
  vcol->SetForceSolid(true);
  lcol->SetVisAttributes(vcol);

  //******************************************************//

  //*****************LP Cylinder**************************//

  //  G4double pulThick        = 0.1*mm;
  G4double pulHeight       = 3.0*mm;
  //  G4double pulOuterRadius  = 4.0*mm;
  G4double pulOuterRadius  = 3.0*mm;
  G4double pulInnerRadius  = 2.8*mm;
  //  G4double pulInnerRadius  = 3.0*mm;
  //  G4double puldistance     = MPMTRadius-1.6*mm;
  G4double puldistance     = MPMTRadius-0.5*pulThick;

  G4Tubs* pul_top = new G4Tubs("pul_top", pulInnerRadius, pulOuterRadius,
			       0.5*pulThick, 0.*deg, 360.*deg);
  G4Tubs* pul_mid = new G4Tubs("pul_mid", pulOuterRadius-pulThick, pulOuterRadius,
			       0.5*pulHeight, 0.*deg, 360.*deg);
  G4Tubs* pul_bot = new G4Tubs("pul_bot", 0, pulOuterRadius,
			       0.5*pulThick, 0.*deg, 360.*deg);

  G4UnionSolid* pul_tam = new G4UnionSolid("pul_tam",pul_mid,pul_top,0,
                                           G4ThreeVector(0,0,0.5*pulHeight+0.5*pulThick));
  G4UnionSolid* pul_tmb = new G4UnionSolid("pul_tmb",pul_tam,pul_bot,0,
                                           G4ThreeVector(0,0,-0.5*pulHeight-0.5*pulThick));

  lpul = new G4LogicalVolume(pul_bot, LP_mat, "lpul");
  //  lpul = new G4LogicalVolume(pul_tmb, LP_mat, "lpul");
  //  ppul = new G4PVPlacement(rotPP,G4ThreeVector(-puldistance,0.,0.), lpul, "ppul", lMFid, false, 0);
  //  ppul2 = new G4PVPlacement(rotPP,G4ThreeVector(puldistance,0.,0.), lpul, "ppul2", lMFid, false, 0);

  G4VisAttributes* vpul = new G4VisAttributes(red);
  vpul->SetVisibility(true);
  vpul->SetForceSolid(true);
  lpul->SetVisAttributes(vpul);

  //******************************************************//  

  return world_phys;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMinEkine
void DMXDetectorConstruction::SetRoomEnergyCut(G4double val)
{
  // set minimum charged particle energy cut - NB: for ROOM
  theRoomMinEkine = val;
  if (theUserLimitsForRoom != 0) 
    {
      theUserLimitsForRoom->SetUserMinEkine(val); 
      G4cout << " Changing Room energy cut to: " << G4BestUnit(val,"Energy")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMinEkine
void DMXDetectorConstruction::SetEnergyCut(G4double val)
{
  // set minimum charged particle energy cut - NB: for Xenon Detector
  theMinEkine = val;
  if (theUserLimitsForDetector != 0) 
    {
      theUserLimitsForDetector->SetUserMinEkine(val);
      G4cout << "Changing Detector energy cut to: " << G4BestUnit(val,"Energy")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMaxTime
void DMXDetectorConstruction::SetRoomTimeCut(G4double val)
{
  // set room time cut:
  theRoomTimeCut = val;
  if (theUserLimitsForRoom != 0) 
    {
      theUserLimitsForRoom->SetUserMaxTime(val);
      G4cout << " Changing Room Time cut to: " << G4BestUnit(val,"Time")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMaxTime
void DMXDetectorConstruction::SetTimeCut(G4double val)
{
  // set detector time cut:
  theMaxTimeCuts = val;
  if (theUserLimitsForDetector != 0) 
    {
      theUserLimitsForDetector->SetUserMaxTime(val);
      G4cout << " Changing Detector Time cut to: " << G4BestUnit(val,"Time")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

//void DMXDetectorConstruction::UpdateGeometry()
//{
//  G4RunManager::GetRunManager()->DefineWorldVolume(Construct());
//}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....


